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ABSTRACT 


This thesis considers the use of acoustic communications in reducing position uncertainty 
for collaborating autonomous underwater vehicles. The foundation of the work relies on 
statistical techniques for accurate navigation without access to GPS, known as 
Simultaneous Localization and Mapping (SLAM). Multiple AUVs permit increased 
coverage, system redundancy and reduced mission times. Collaboration through acoustic 
communications can minimize navigational uncertainty by permitting the group to benefit 
from locally discovered information. However, the propagation of acoustic 
communications can be used to counter detect the system during naval operations. 

The thesis gives explicit consideration to tactical security in acoustic 
communications for a multi-AUV SLAM system. It provides initial techniques and 
analysis for minimizing communications between AUVs. The reduction is accomplished 
through a statistical method that allows for the estimation of the updated covariance 
matrices. Normally, SLAM techniques use expropioceptive (sonar and cameras) sensors 
and computer vision algorithms for the detection and tracking of navigational references. 
We propose a novel use of the acoustic modem as another sensor. It leverages the 
physical characteristics of underwater acoustic transmissions and the information 
transmitted in the signal to provide an additional measurement. We believe this is the first 
emphasis on minimizing communications within a multi-vehicle SLAM approach. 
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I. 


INTRODUCTION 


A. MOTIVATION 

Few naval powers possess the ability to challenge the United States Navy in deep 
water or effectively deny us access to anywhere we wish to sail. The rise of cheap, 
effective, asymmetric anti-access and area denial (A2/AD) systems represents the best 
possible solution to denying the United States access to a region without the need to be 
able to conventionally confront maritime forces in Mahanian naval combat. The United 
States experienced this problem first hand on April 14, 1988, when the USS Samuel B. 
Roberts (FFG 58) struck a mine in the Persian Gulf while escorting tankers as part of 
Operation Earnest Will. The subsequent transport back to the United States and repairs 
cost taxpayers $89.5 million [1]. 

The mere threat of a naval minefield can effectively stop maritime traffic and 
commerce in a port or strategic chokepoint, such as the Strait of Hormuz where the USS 
Samuel B. Roberts was operating. The low cost and advancing technology of naval mines 
makes them particularly well-suited for A2/AD applications. The United States must 
possess the ability to rapidly and covertly map and neutralize a minefield in order to 
assure access and ensure the free flow of maritime trade in the global commons. The 
present means of mine countermeasures largely reside on surface ships and are neither 
covert nor rapid. Autonomous underwater vehicles (AUV) have the ability to map the 
environment in a covert manner that does not necessarily require an overt presence. They 
do not require real time human control and are difficult to detect. Multiple AUVs permit 
a wider area to be searched and mapped more efficiently than any present means. The use 
of multiple AUVs will be required to keep minefield mapping and clearance both covert 
and efficient. At present, the use of AUVs in minefield mapping is in its infancy and 
AUVs lack the ability to collaboratively complete a mission; each vehicle would operate 
essentially independently. Recent research has enabled multiple AUVs to begin 
coordinating their efforts. This body of research has significant operational value for 
undersea operations. 
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B. PROBLEM STATEMENT 


The underwater environment presents significant challenges to both navigation 
and localization since the vehicles must operate without the benefit of the global 
positioning system (GPS). This environment necessitates a more accurate means of 
navigation for AUVs to ensure that they can remain underwater and undetected. 
Obtaining a GPS fix for an AUV is a highly inefficient part of the overall mission profile 
as it necessitates the AUV rising from deep water, loitering on the surface where the 
threat of collision or counter-detection is greater, submerging back to its programmed 
search depth, and reacquiring its position with respect to the underwater environment. 
Underwater beacon systems have been developed and fielded, but deployment of these 
beacon fields requires both prior warning and time that may not be available. AUVs must 
have the ability to localize their position without the need for external navigation aids. 
The development of simultaneous localization and mapping (SLAM) [2] allowed 
autonomous vehicles to both localize their position and map their environment at the 
same time. This field provides a major operational capability for accurate underwater 
navigation and mapping. Many of the Navy’s undersea warfare missions would benefit 
from advances in SLAM, including mine countermeasures. Research into SLAM 
algorithms has exploded in recent years as researchers and end users search for a way to 
make these robots reliant on their onboard systems only and make them more capable to 
accurately navigate in difficult terrain and environments. However, only a small 
percentage of the field has tried to solve the problems of multiple AUV coordination in a 
SLAM environment. Collaborative SLAM between multiple AUVs permits improved 
coverage, greater accuracy, and faster, more efficient operations. In the given minefield 
mapping scenario under current employment constructs, a fleet of AUVs tasked to map a 
minefield would be operating in the same space, but independently. The reduction of 
their position uncertainty from either SLAM or a GPS fix cannot be shared with other 
vehicles for their benefit. The ability for a single AUV to share its position, or state, 
information following a reduction in its position uncertainty, whether through a GPS fix 
or through SLAM, with the other AUVs conducting SLAM operations would prolong 
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submerged operations and minimize disruptions to mapping operations while potentially 
improving the overall quality of the maps produeed. 

This thesis will eonsider the role of aeoustie eommunieations in multiple-AUV 
operations. Multiple-AUV SLAM (MVSLAM) is the ability for multiple AUVs to 
simultaneously eonduct SLAM operations and share their information with other AUVs 
to improve the overall performanee of the group. This thesis will explore the value of 
aeoustie eommunieations to MVSLAM in reducing position uncertainty, and tangentially, 
map accuracy. Since the acoustic communications present a real threat to counter¬ 
detection, this thesis will heuristically balance performance improvements with tactical 
security considerations. It will explore strategies for minimizing communications while 
maintaining navigational accuracy. 

Tactical security concerns, an issue unique to the military, will underpin this 
thesis to ensure that the developed solution keeps with the strong desire to remain covert 
and undetected. This viewpoint is conspicuously absent from the present body of research 
and will be a major contribution to it. Current approaches rely on frequent 
communication between vehicles to pass information; whether it is position information, 
maps, or command and control functions is immaterial. The frequencies of the acoustic 
modems used in this thesis will propagate, under typical sound conditions, omni¬ 
directionally for several kilometers and correspond to frequencies that active sonar 
intercept receivers can detect. Minimizing the number of messages required to be 
transmitted and the intervals at which they need to be transmitted directly correlate with a 
reduction in the probability of counter-detection. 

C. APPROACH 

The overall objective of this thesis is to adapt the Second Generation Incremental 
Smoothing and Mapping (iSAM2) algorithm [3] for use in multiple-AUV SLAM. To 
accomplish this, this thesis will be organized as follows. Chapter II will describe the 
equipment used in this research and review the literature of sound propagation in 
seawater and acoustic ranging and localization methods. Chapter III will review current 
approaches to MVSLAM, the fundamental mathematics behind position uncertainty in 
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robotics, and provide a detailed deseription of iSAM2. Chapter IV will discuss the 
development of a hybrid approach to MVSLAM using the mathematieal concept of a 
Bayesian inferenee eoupled with a novel use of aeoustie eommunieations to reduee 
position uneertainty. We will show that this solution ean be applied to n-number of 
AUVs. Chapter V will validate these ehanges through simulation, along with an 
exploration of the value of aeoustie eommunieations to MVSLAM. Taetieal seeurity 
eonsiderations will be heuristieally ineorporated. Chapter VI will diseuss the major 
results, eontributions, and eonelusions from this work, as well as propose signifieant 
areas for future work. 
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II. ACOUSTIC COMMUNICATIONS AND RANGING 


Communication is fundamental to the suecess of MVSLAM operations. In the 
underwater environment, only acoustie eommunieations have the necessary range to 
enable eommunieations between vehicles. Underwater eommunication is challenging 
beeause of the highly uneertain nature of sound propagation, the low data rate in the 
eommunieations ehannel, and the signifieant losses the signal ineurs during interaetions 
with water moleeules. Despite these ehallenges, understanding how aeoustie rays 
propagate through the water provides information that ean be used in improving 
navigational aeeuraey. This seetion will discuss the equipment used in this thesis, the 
propagation of sound in the oeean, and a method of acoustie ranging. 

A. MODIFIED REMUS-IOO AUV 

This thesis exelusively utilizes the Hydroid REMUS-100 AUV for SLAM 
operations. NFS owns two REMUS AUVs and has modified them extensively to support 
various researeh aims. This seetion will detail both the general eharaeteristies of the 
REMUS 100 AUV as well as the speeifie modifieations that NFS has made to the 
vehieles. Figure 1 shows the NFS REMUS AUV during a photo opportunity during the 
National Aeronautics and Space Administration (NASA) Extreme Environments Mission 
Operations (NEEMO) off the eoast of Key Largo, Florida, in September 2013. The 
mission used the REMUS vehieles to map a simulated asteroid environment and a 
tethered hovering AUV as a guidanee and astronaut-assistance platform. 
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Figure 1. Two NASA astronauts hold the NPS-modified REMUS 100 AUV 
during a photo opportunity at the Aquarius Reef Base during the NEEMO 
mission off the coast of Key Largo, Florida, in September 2013. 


The REMUS 100 is a man-portable, lightweight AUV designed primarily for 
survey-type operations. Table 1 outlines the basic specifications. 


Table 1. Specifications and operating characteristics of the NPS-modified 
Hydroid REMUS-100 AUV with the BlueView forward looking 
sonar and cross body thrusters attached. 

Specifications of the REMUS-100 AUV 


Length 

Approx. 8 ft 

Diameter 

7.5 in 

Weight 

95 lbs 

Maximum Depth 

100 m 

Speed 

Up to 4.5 kts 

Endurance 

8-10 hours 
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The REMUS 100 is a highly modular system that can be easily customized to 
complete a wide variety of underwater tasks. For the type of research that NFS conducts 
with these vehicles, we have the following sensors and systems installed: 

• Kearfott KN-6051 SEADeViL inertial navigation system (INS) with a 
Doppler Velocity Eog (DVL) and GPS 

• Fore and aft cross-body thrusters from Hydroid 

• Woods Hole Oceanographic Institute (WHOI) Acoustic Micromodem 

• YSl CT-600-XE Conductivity-Temperature-Depth (CTD) Sensor 

• Marine Sonic Technology, Etd. Side Scan Sonar 

• Teledyne BlueView MB2250 3D Microbathymetry Sonar 

• Teledyne BlueView FE450X 2D Forward Fooking Sonar (FES) 

• Acoustic Doppler Current Profiler (ADCP) 

The remainder of this section will describe the INS, acoustic modem, and forward 
looking sonar in greater detail since they are fundamental to this thesis. 

1. Kearfott KN-6051 SEADeViL INS 

The KN-6051 is a military-grade INS that combines inputs from onboard sensors, 
such as the DVF, with external measurements from GPS when they can be acquired. The 
ring laser gyro-based system has an overall accuracy of 0.5% error per unit distance 
traveled. This equates to a 5-meter error per kilometer. This error is in terms of the 
circular error probable (CEP) rate, which translates to a 50% probability of being within 
that circle. For heading, the INS is accurate to 5 mils, or 0.28 degrees. The GPS, with 
DVF aiding, is only accurate to a 10-meter CEP [4]. However, since it is an external 
position source, it is a community standard practice to assume no error in that position 
and take it as truth. 

2. WHOI Acoustic Micromodem 

WHOI developed this acoustic micromodem for their research purposes and it has 
quickly become an industry-leading piece of equipment. The device requires very little 
power. It idles and receives at approximately 158 mW on the 12V system, and requires 

less than 100 W for a five second burst transmission—very low power. The modem 
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transmits at approximately 25 kHz. Acoustic communications in seawater are very low 
data rate channels compared to electromagnetic or radio communications in air, and are 
on the order of bytes to a few kilobytes per second [5]. 

3. Teledyne BlueView FL450X FLS 

The sonar systems installed come as a specially-engineered module from 
Teledyne BlueView to include both the FLS and microbathymetry sonar in the same 
housing. Since this thesis does not require the use of microbathymetry information, we 
will discuss only the performance of the FL450X FLS. The sonar operates at 450 kHz 
and has a field of view of 130 degrees horizontally and 45 degrees vertically. The sonar 
can detect objects in this field of view out to 280 meters, but between 5-100 meters is 
optimal. The object detection software supplied with the sonar can detect objects with an 
accuracy of 1 m in the range direction and 1.2 degree accuracy in bearing [6]. 

B. SOUND PROPAGATION IN THE OCEAN 

1. Overview 

The ocean is an incredibly complex environment and sound propagation in the 
ocean can be exceedingly difficult to model and predict. Simplifications and assumptions 
about sound propagation in one area of the world may not hold true for another area 
simply based on physical conditions such as bottom type, presence of biological 
organisms, salinity content, and more. However, basic models of sound propagation and 
sound speed equations will function adequately over the limited physical ranges that 
multiple-AUV SLAM will encompass. This section will discuss the basic mechanisms of 
sound propagation, the factors that influence sound speed, and the factors affecting 
transmission loss from a radiated source. 

2. Sound Speed 

Three principle factors govern sound speed in water: salinity, pressure, and 
temperature. In the littoral waters where minefields would likely be placed, salinity and 
temperature dominate the sound speed profile. Pressure does not normally affect the 
change in sound speed until the depth is below the main thermocline, which generally 
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occurs only in deep water. Several researehers have attempted to ereate an equation to 
prediet the speed of sound in seawater. The equation developed by Del Grosso [7] and 
updated by Dushaw et al. [8]7 has been aeeepted as the most accurate and useable 
equation by the aeoustics eommunity. Equation (2.1) shows the Del Grosso equation with 
the modifieations from [8]. The numerieal eonstants (Cxx) for eaeh term are provided in 
Appendix A. The physieal properties are measured in degrees Celsius for temperature 
(T), parts per thousand (ppt) for salinity (S), and kg/em^ (gauge) for pressure (P). 

C^jp\iTi / 5] = Cqqq + ACj. + AC^ + ACp + AC^pp 

1402.392 

ACp = Cp^T + Cp^ + Cp-^T 

ACs=CsiS + Cs2S' (2.1) 

ACp = Cp^P + CppP + Cp^P 

AC^PP = CpgTS + CppTP + Cp2p^ P '^Cpp^P + Cpp^TP +... 

Cp2pT P + Cg2P2^ ^ ^^TSIP^^ ^ + CpgpTSP 

The equation performs extremely well in both deep and shallow water. Reported 
aceuracy, verified in [8], shows the equation to be aeeurate to within 0.3 meters per 
seeond (m/s) with a standard deviation of 0.05 m/s aeross the range of likely input values 
for temperature, salinity, and pressure. Equation (2.1) will be used onboard the REMUS 
vehicle to calculate sound speed as sensor measurements become available for the 
purpose of aeoustie ranging. 

3. Sound Propagation in the Ocean 

In an isoveloeity sound profile, sound propagates linearly and spreads in a 
spherieal manner until it internets with a boundary layer, sueh as the air-oeean or oeean- 
bottom interfaees. However, isoveloeity sound profiles are unlikely to exist in a dynamie 
ocean environment; and therefore, the sound propagates in a curvilinear manner, 
refracting incrementally as governed by Snell’s Eaw, given below in Equation (2.2), 
where c is the sound speed in a given layer and 0 is ineident angle of the sound ray [9], 
[10]. The eurvilinear behavior results from ineremental ehanges in temperature, pressure. 
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and salinity in the vertical water column, as a qualitative examination of Equation (2.1) 
shows. 


008 6*; _ 008 6*2 
Cl Cj 


( 2 . 2 ) 



Figure 2. Depiction of Snell’s Law where Cj and Cj are the sound speeds in the 
given water layer and 0^ and 0^ are the grazing angles of the acoustic rays 
from the horizontal plane dividing the two water layers in question. 

Figure 2 and Equation (2.2) show the relationship between how changes in the 
speed of sound in the vertical water column produce curvilinear ray propagation. Figure 3 
shows how sound propagates horizontally for several non-isovelocity sound speed 
profdes. 
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Figure 3. A collection of sound propagation plan views for varying sound speed 

profdes, from [11], 

Over the ranges that AUVs with the WHOI acoustic modem can communicate, 
the propagation paths will be either curvilinear direct path, with no boundary interactions, 
as would likely be the case in deep water, or will reflect off the bottom, surface, or both, 
resulting in greater losses and shorter ranges, as we would expect in a very shallow water 
littoral environment. Snell’s Law provides a convenient way of thinking about sound 
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propagation at the conceptual level, but the formula does not capture sufficiently the 
complexity of the propagation. The differential equation form, independent of frequency, 
can be iteratively solved through computer-based numerical methods, such as finite 
element analysis. 



Scund speed (m/s) 



.1101 -^^- 1 - 1 - 1 - 1 

0 0.5 1 1,5 2 2.5 3 
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Figure 4. On the left (a), the canonical Munk sound speed profile, ending at 100 
meters depth. On the right (b), the ray trace diagram from a source radiating 
omni-directionally in (a) at 30 meters depth with no accounting for 
transmission losses. The red and blue rays simply indicate whether the initial 
transmission angle was above or below the horizontal plane [12]. 


Figure 4 shows how busy the acoustic picture can become with an 
omnidirectional source radiating in the ocean at a given depth. Each ray shown is an 
individual ray transmitted from the source at a specified angle. The sum of these rays 
reflects the behavior of an omnidirectional source. The ray traces provide no information 
regarding the strength of the acoustic signal at any point since transmission loss varies 
with multiple parameters. Since the degradation of an acoustic signal and the length of 
time it takes to propagate from the source to the receiver are of vital interest to this thesis, 
a discussion of the sonar equation and its subordinate loss terms is therefore warranted. 
Additionally, Figure 4 provides an initial insight into the uncertainty that we face in 
developing acoustic ranging equations, given the inability to know which ray we 
received, how that particular ray interacted with the ocean bottom and surface, or how the 
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local effects of the seawater it passed through affeeted it, to name a seant few of the 
variables at play here. 

4. The Passive Sonar Equation 

In its simplest form, the passive sonar equation relates the received signal level to 
the level required for signal deteetion and is generally reported in deeibels (dB), 
refereneed to 1 ^Pa in seawater. i From source to operator, the signal undergoes a 
number of losses. For ease of representation, we present them in a tabular format in Table 

2 . 


Table 2. The major terms of the sonar equation and their definitions [9], [10]. 


Parameter Symbol 

Description 

SL 

Source level: the radiated intensity of 

the acoustic source, referenced to 1 

yard from the source, by convention. 

TL 

Transmission losses: includes 

spreading and attenuation, one-way 

NL 

Noise level: includes ambient and self¬ 
noise at the receiver 

DI 

Directivity index: the ability of the 

array to detect the signal 

DT 

Detection threshold: the ability of the 

system to detect the signal 


^ All signal values in the remainder of this paper will be referenced to 1 JLlPa in seawater, so we shall 
drop the additional text. 
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These terms eombine to form the passive sonar equation [10]: 

SL-TL>NL-DI + DT (2.3) 

The right hand side of Equation (2.3) relates the teehnieal performanee of the 
sonar system in its immediate environment, in this ease the acoustic modem receiver on 
the AUV in seawater. The detection threshold and directivity index, closely tied to signal 
processing, are immaterial to this thesis and will not be discussed further. The noise level 
consists of the noise at the receiver, which comes from two parts: first, the self-noise 
within the receiver itself, and, second, the ambient or background noise at the receiver 
from biologies, shipping, seismic movements, weather, etc. At the frequency of the 
WHOl acoustic modem, 23-27 kHz, the background noise is primarily wind-driven with 
some biologic activities [13]. The source level of the WHOl acoustic modem is a fixed 
quantity related to the lOOW transmission power, which equates to approximately 190 dB 
[5]. The remaining term, transmission losses, consist of spreading and absorption effects 
and losses at the boundaries. They will occupy the remainder of this section since they 
drive the operating range of the acoustic modem and form the core of acoustic ranging 
methods. 

5. Spreading and Absorption 

Sound radiates from the acoustic modem in a roughly spherical manner until it 
interacts with a boundary layer, whether the surface or the bottom. While travelling, each 
of the rays follows a curvilinear path as previously discussed in Section 11.B.3. If, at the 
point of transmission, we consider the acoustic message to be a sphere with a finite 
amount of energy, the energy at the wave front decreases proportionally to the increase in 
the surface area of the sphere, which can be approximated as the square of range (r). The 
same amount of energy must cover this larger area and thus the intensity weakens. 
However, if we assume that we are operating in shallow water and that the acoustic 
transmission will interact with a boundary layer, the spreading becomes cylindrical and 
reduces the loss rate, in dB, by half, as shown below [10]. 
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TLcyiindrical = 1 0 lOg T 


(2.4) 


Absorption results from three separate phenomena: shear viseosity, ionie effeets, 
and pressure (depth). Transmission loss due to absorption is substantially worse in 
seawater than in pure water due to the ionie equilibria. We need not dive into eaeh 
phenomenon in detail, but the cumulative effects of absorption at the frequency of 
interest amount to 5-10 dB per kilometer, assuming a constant pressure (depth) and 
depending on the temperature and salinity of the seawater [10]. 


6. Losses at the Boundary Layer 

As Figure 4 showed, a sound signal from the acoustic modem can reflect off the 
surface and bottom multiple times before reaching a receiver. Each of these bounces will 
incur losses. At the air-ocean interface, nearly all the sound is reflected, vice transmitted 
through the boundary, so the loss is smaller, though rough surfaces can result in 
additional losses from scattering. At the ocean bottom, the losses are governed primarily 
by Snell’s Law, including the ray grazing angle {0 ), sound speed (c), and the density 
(p) of both the ocean floor and water layer directly above it [14]. 

BL = 2Q\og,,\R,,\ 

_ msm0^-{n -cos^ 

“ m sin ^i+(n'-cos'll)'"' 

Pi ^1 
m = —^,n = — 

P\ 

The porosity of the ocean bottom and the grazing angle drive the losses, with 
additional input from the sedimentary layers, surface roughness, presence of biological 
materials, air bubbles, and other attenuation effects. Providing a single, shorthand 
calculation, as we did for cylindrical spreading losses, is not possible. However, Urick 
[10] reports that bottom losses can range from 0-30 dB per bounce. 
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C. ACOUSTIC RANGING, ONE WAY TRAVEL TIME 

One-way travel time (OWTT) is the most widely aecepted model for acoustic 
ranging in the field today [15]-[19]. This method reduces the complexities presented in 
the preceding sections into the standard kinematic formula of speed multiplied by time. 
The advances in the field have largely revolved around the timing aspect of the equation, 
through characterization of precision clocks, as in [17], and the programmatic and 
algorithmic details of reducing uncertainty in the actual time of flight (TOP) 
measurement. Discussions of the sound speed aspect of the equation have remained 
largely unaddressed, with researchers assuming an isovelocity sound speed profde and, 
therefore, straight line, vice curvilinear, ray paths. These ranges underestimate the actual 
range because of the effects of Snell’s Law as previously discussed. However, the 
accuracy is within acceptable margins (<1 m) for the range of interest and the 
calculations are not computationally expensive to run, so the OWTT framework makes 
logical sense. Accounting for transmission losses in acoustic ranging explicitly adds 
several orders of magnitude in difficulty given the variability of the factors to weather, 
geography, salinity, bottom type, etc., and therefore will not be addressed in the OWTT 
formula or further in this thesis. 

The OWTT equation will be calculated using the standard kinematic formula: 

Range = ) (2.6) 

Calculating the time of flight will be accomplished through the timestamps that 
the WHOI acoustic modem applies to the incoming and outgoing messages. The speed of 
sound will be calculated from the corrected Del Grosso sound speed equation discussed 
in Section II.B.2 using measurements from the onboard CTD sensor. In terms of 
uncertainty calculations, empirical testing of the clocks onboard the REMUS vehicles at 
NPS indicates no clock drift over eight hours. Similar results from [16], [17], [19], and 
[20] validate this conclusion. Therefore, we will neglect any uncertainty added from time 
synchronization issues. For the uncertainty of sound speed, the results from [21] that 
show the Del Grosso sound speed equation to have an assumed Gaussian distribution 
with a mean of 0.3 m/s and a standard deviation of 0.05 m/s, assuming perfect sensor 
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inputs for pressure, temperature, and salinity. To put the uneertainty into physieal 
perspeetive, two AUVs 500 meters apart would see an aeoustie range error of 0.0996 ± 
0.0332 meters to two standard deviations (2cr). Integrating this equation into the system 
model needed for MVSLAM will be done in Section IV. 

This section presented the equipment used in this thesis, reviewed sound 
propagation in an ocean environment, and discussed a simple, widely-used technique for 
acoustic ranging. We can use the highly uncertain nature of sound propagation to help 
reduce position uncertainty in an AUV. The implementation of that concept is the focus 
of both the existing SLAM literature and algorithms, presented next in Section III, and in 
the proposed use of acoustic communications for navigational accuracy developed later 
in Section IV. 
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III. SIMULTANEOUS LOCALIZATION AND MAPPING 


This section builds the case for probabilistic robotics from first principles by 
examining adaptation of stochastic elements into traditional control system formulations. 
It will then discuss, in depth, the development of iSAM2, the SLAM algorithm used on 
the REMUS vehicles at NFS, and will illustrate the performance of the algorithm with a 
small example. iSAM2 will play a foundational role in the development of a distributed, 
collaborative SLAM framework in Section IV. 

A. INTRODUCTION 

1. Position Uncertainty in Robotics 

Smith, Self, and Cheeseman [22] first postulated positioning in robotics as a 
stochastic, rather than deterministic, problem. This fundamental shift in perspective 
allowed engineers and scientists to incorporate positional uncertainty in the design of the 
robotic system instead of simply coping with the degradation in quality that could result. 
The incorporation of probability into this process stems from the realization that no 
sensor or system is perfect in its sensing or movements, respectively. For example, with 
the REMUS 100 vehicle, the BlueView forward looking sonar can only sense objects to 
an accuracy of 0.1 meters in the range direction and 1.2 degrees in the bearing direction. 
The onboard control system can command the electrical motor to spin the propeller at 
2500 RPM, but the known and unknown characteristics of the system may result in the 
propeller only spinning 2480 RPM, for example. The difference results in slightly less 
speed, which forces the uncertainty about the vehicle’s position to grow. The 
mathematical representation of a system expanded from the dynamics model to include a 
system covariance matrix that shows the probabilistic relationship between the state 
variables. 

To ground the work, we discuss the approach to SLAM algorithms from its roots 
in control systems theory. In Equation (3.0), shown in the canonical discrete state-space 
form, the first row constitutes the system process or dynamics model, and the second the 
measurements by the system [23]. These are matrix equations resulting in assumed linear 
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relationships. The position at the next time step, is a funetion of the dynamieal 
model, sueh as dead reekoning navigation of an AUV, with additive noise co^. Unless 

explieitly known to be otherwise, all system noise is assumed to be independent and 
identieally distributed (i.i.d.), zero-mean and Gaussian, also known as Gaussian white 
noise. The measurements of the system, z^, from aeoustie ranging via a beaeon, for 
example, are also subjeet to Gaussian white noise v^. An illustration of this in the 
eontext of SLAM will be benefieial. 




(3.0) 


2. Sources of Uncertainty 

The eombination and propagation of the proeess and measurement noise, and 

respeetively, results in a probabilistie distribution of the vehiele’s position. That 

distribution is assumed to be Gaussian sinee all eomponents of the uneertainty terms are 
Gaussian as well. The vehiele has a ealeulable probability of being at a given point from 
the mean position. The same ean be applied to the position uneertainty of landmarks. For 
ease, we define the position uneertainty (PUC) as a Gaussian ellipsoid and it is the 
graphieal estimate of the position uneertainty of the vehiele at a given time. The ellipsoid 
ean be projeeted onto the horizontal plane, as an ellipse, to show the area that the AUV 
may be in. We assume that all ellipses represent a 95% probability (two standard 
deviations) of being inside the ellipse. The PUC size ehanges over time as the errors 
eompound. 

From an inertial navigation-only perspeetive, the PUC will grow over time as the 
errors in the INS eompound. These errors stems from the engineering toleranees of the 
gyros and aeeelerometers in the INS, the ability to preeisely sense the position of the 
inertial measurement units or motion of the aeeelerometers, and the preeision with whieh 
the software ean extrapolate the estimated position from that noisy data. The inputs into 
that ealeulation, sueh as the eommanded vehiele speed, gyroeompass heading, the eurrent 
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data from the ADCP, etc., all introduce additional errors. The magnitude of the error for a 
given time step is relatively constant, usually varying with physical parameters, such as 
current, but over time, those errors are additive. For the Kearfott INS installed on the 
REMUS, the error can be summarized as approximately 0.5% per distance traveled, as 
previously stated. In Equation (3.0), these uncertainties manifest in the 64 term. 

Since we will be examining SLAM, which relies on the ability of the vehicle to 
sense environmental features, the errors within the sensing systems must also be 
considered. Eor example, a forward looking sonar mounted on the front of the AUV 
provides image data up to 5 Hz relative to the AUV’s position. A detection and tracking 
procedure can be used to resolve the relative bearing and range to the vehicle’s nose. That 
procedure will have uncertainty associated with it, manifesting as uncertainty in bearing 
and range, which can be transformed to the AUV’s local coordinate system for 
incorporation. The errors from the sensing equipment reside in the term of Equation 
(3.0). A small example in the next section provides additional insight into this area. 

3. SLAM as a Stochastic Process 

SLAM, as the name suggests, is the ability to map an area and navigate off of that 
map in real time. With the implementation of SLAM onto a mobile robot, we can 
consider the robot conducting two primary evolutions to accomplish this: sensing the 
environment and movement within that environment. In some robotic applications, the 
robot must stop to sense, while in others, such as underwater, the vehicle is continually in 
motion while sensing, to include hovering or station-keeping. However, the discretization 
of the latter case removes the complications associated with that motion. 

Using the REMUS vehicle as an example, we can see how these probabilistic 
relationships function in a SLAM environment. The REMUS vehicle uses a twelve 
variable state vector, or pose (x), in Equation (3.1) that provides six degrees of freedom 
in movement since it contains both linear (x,y,z) and angular {(f), 9,y/) position and rates 
(u,v,w for linear, and p,q,r for angular). The associated covariance matrix, or relationship 
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between each of the state variables, is, by definition, a 12x12 matrix. Or more generally, 
an n-dimensional state vector will have an n-by-n covariance matrix, E(x). 
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The off diagonal terms in the covariance matrix represent the cross-correlation 
between the individual state variables. 

The robot begins by sensing a landmark from its current position. At this point, 
the robot has defined a PUC, as indicated by the red circle in Figure 5. 



Figure 5. The robot, with a finite PUC, senses a landmark and takes a sensor 
measurement. The uncertainty of the feature and the AUV are the blue ellipse 

and red circle, respectively. 
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The sensor measurement of the landmark will have uncertainty associated with it, 
as previously stated. In the simplest sense, the error from the measurement of the 
landmark can be superposed onto the vehicle and combined with the vehicle’s current 
PUC through a Bayesian inference to produce a new PUC. The mathematical 
implementation of the Bayesian inference will be covered in Section IV.C.l. The vehicle 
then moves and takes another measurement of the same landmark, shown in Figure 6. 





Figure 6. The robot, with an updated PUC, has moved a fixed distance and takes 
a measurement of the same landmark. The initial PUC is in light red for 

comparison. 


The incorporation of the sensor measurement is again superposed onto the vehicle 
and algorithmically combined to form a new PUC measurement. In Figure 7, we see the 
result after the second time step. 
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Figure 7. After two measurements of the same landmark, the robot’s PUC, the 
red ellipse is substantially smaller than it was at the beginning. The initial 
ellipse is in light red to show the reduetion. The updated uneertainty of the 
feature is shown as the light blue ellipse, and is markedly smaller than the 

initial, shown in dark blue. 

We can see visually that using sensor measurements of environmental features 
allowed us to reduce the AUV’s PUC significantly. Simultaneously, all measurements of 
the feature are integrated through a Bayesian inference or other statistical methods. The 
resulting positional uncertainty of the AUV and the feature is smaller than if only a single 
measurement had been taken. The end result of the SLAM process is that the AUV has 
constrained its navigational uncertainty, which improves map accuracy, while creating a 
map of the environment. 

This simple example showed the basic process for SLAM. The effects on each of 
the twelve individual state variables may or may not be related to the SLAM process, 
depending on the information obtained with the measurements. In most implementations 
in the underwater environment, the SLAM process only affects linear position in the 
horizontal plane (x,y) and heading (i// ). The rest of the state variables are generally a 
function of the vehicle control systems and physical system design. For notational 
convenience, they will be omitted going forward in this thesis. 

This example also highlights the challenges associated with SLAM. While not 
discussed in the context of the example, several other algorithms are vitally important for 
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SLAM operations. First, the deteetion and traeking of features is not an exact science, as 
discussed in Section III.A.2. Increasing the performance of that algorithm will improve 
navigational and map accuracy. Second, this example assumed the ability to 
unambiguously differentiate between features. Only one feature was presented in this 
example, but features may appear close together and associating measurements with the 
correct feature is crucial to ensuring stable algorithm performance and maintaining map 
accuracy. This process of correlating measurements with features is known as data 
association. Lastly, the ability to return to a previously detected feature and recognize it 
as previously detected is known as loop closure. Loop closure is a vital aspect of making 
SLAM algorithms with operational utility. 

4. Optimal Estimation 

Rudolf Kalman’s [24] seminal paper in 1960 proving the optimality of his new 
linear filter opened a new field known as optimal estimation, in which we attempt to 
estimate uncertain processes such as AUV motion. Kalman filters assume the system is 
linear and the noise is Gaussian i.i.d and optimally balance the dynamic uncertainty 
inherent in process and measurement models. The linear Kalman filter has spun off into 
several new fields of filtering: particle filters, extended Kalman filters (EKF), which 
considers non-linear motion, extended information filters, and smoothing filters, to name 
a few of the major ones. Most robotic systems use an EKF, or a variant thereof, in SLAM 
for computational efficiency. But some, such as the filter implemented on the NFS 
REMUS vehicles, use a smoothing filter. 

EKEs consider the information at a given time step, with all previous information 
summarized in the prior state estimate. The EKE incorporates any sensor measurements 
and linearizes the equations of motion at the current time step to produce a updated 
estimate of the robot pose. Smoothing filters, on the other hand, consider all prior 
information at each time step without summarizing. While they have historically been too 
computationally inefficient for most real time applications, faster smoothing filters have 
emerged concurrently with processing power and hardware advances, iSAM2 being one 
of them. Using iSAM2 provides a more accurate estimate than a traditional EKE-based 
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SLAM approach since it revisits all previous estimates of features and navigation 
estimates at eaeh time step, which is highly desirable in the underwater environment. 


C. INCREMENTAL SMOOTHING AND MAPPING 


1. Process and Measurement Models 

Developed by Kaess and Dellaert [3], [25], iSAM2 fulfdls the eore goals of any 
SLAM proeess: exaetness and eomputational effieiency. An EKF, essentially a non-linear 
Kalman Filter, linearizes about a point at a given time step, whieh is then ineorporated 
into the prior state estimate, unable to be ehanged in the future. iSAM2 provides for fluid 
variable relinearization, which seems eomplex, but is very efficient sinee in SLAM 
applieations, the information matrix underlying the feature map that we are trying to 
linearize about is at all times sparse [25].^ 

A SLAM proeess is affeeted by both the executed vehiele trajeetory and a map of 
the features, whieh as previously stated, are both uneertain. Kaess starts with the 
assumption of a non-linear model and eonverts this into a least squares problem, whieh 
allows us to estimate all unknowns given the available measurements through a maximum 
a priori (MAP) estimate. We begin by stating that SLAM follows a basie proeess and 
measurement model, akin to the state-spaee model provided in Equation (3.0): 

x = Ax + Bu + CO 

(3.2) 

z = Hz+ 0 


ryand u are the zero-mean, i.i.d. Gaussian noise assoeiated with proeess and 
measurements, respeetively. In this ease, the A and B matriees, the equations eoneerning 
dynamies and control inputs, can be collapsed into a single function that describes the 
dynamie motion of the vehicle—the proeess model. The measurement equations ean be 
similarly adjusted to aeeount for the measurement of landmarks and the subsequent data 
assoeiation. Equation (3.2) thus reduees as shown below: 

x = f{x,u) + co 
z = h{x,l) + V 

2 All information and equations in this section, for subsections 1-5, are from [25] unless explicitly 
cited otherwise. 
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I is the landmark that has been sensed and faetors into the measurement model, z, 
of Equation (3.2). The explieit time dependeneies of the system have been omitted for 
simplieity, known as diseretization, but these equations do vary with time, and multiple 
landmarks ean exist. 


2. Linearization 


The eonstruetion of the MAP estimate through least squares requires the 
linearization of both the proeess and measurement models around the current estimates. 
We accomplish this through a Taylor series expansion of the two equations in Equation 
(3.3), which yields the resulting three Jacobian matrices. 


F = 


df(x,u) 


dx 


\xf ,4 


dh(x,l) 


J = 


dh(x,l) 


dl 


(3.4) 


These matrices are collected into a larger matrix A along with a fourth special 
matrix, G, that allows us to not consider the terms. These matrices are oriented in 
the A matrix as follows: 


A = 


GOOD 
F H 0 0 
0 F H 0 
0 7 0 7 


(3.4) 


The least squares process can be qualitatively described as seeking the minimum 
of a given argument consisting of squared terms. In this case, the four Jacobian matrices 
and the navigation (process) and measurement prediction error terms, aandc, 
respectively, can be arranged to represent the least squares estimation problem: 


S0 = arg min 

se 


M 

Ilk' 


-1 


Sx.^^ +G\5x- - a- 




k=\ 


(3.5) 
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After some algebraie manipulation, it ean be proven that Equation (3.5) reduees as 
follows: 

=argmin||A6>-Z7f (3.6) 

Equation (3.6) is now in the standard form for linear least squares (EES) 
estimation, where 6 is the eoneatenation of the vehiele pose and the landmark variables 
and b is the eoneatenation of the navigation predietion and measurement predietion 
errors. 9* beeomes the new predietion. This EES problem ean be solved using standard 
methods, sueh as QR faetorization or Cholesky deeomposition. Kaess opts for QR 
faetorization sinee it paves the way for eomputationally easier ineremental updates to the 
A matrix, whieh grows over time. The unique solution to the EES problem is also termed, 
in this applieation, the square root information matrix (SQIM). iSAM2 updates the SQIM 
as new information beeomes available. It gains effieieney by using the previous solution 
and only performing the ealeulations on the new measurements. 

3. Variable Reordering 

Variable reordering is a linear algebra teehnique that reorders the eolumns of the 
information matrix. Kaess applies it in bloeks, eaeh of whieh eonforms to a single pose 
node or landmark [3]. This method reduees the duplieation of landmarks if the robot 
revisits a previously sensed landmark, or eloses the loop in eommunity standard 
language. Tailing to reeognize loop elosure events in SEAM problems leads to 
duplieation of previously visited landmarks, and, for iSAM2 in partieular, the appearanee 
of non-zero entries in the R matrix following QR faetorization—a highly undesirable 
result. Kaess solves this problem through the applieation of a eommon linear algebra 
teehnique, eolumn approximate minimum degree (COEAMD), to bloeks in the R matrix 
that eorrelate to previous AUV positions and landmark positions. 

4. Process Results 

The proeess deseribed in the preeeding seetions is eondueted eaeh time the AEIV 
senses a new landmark. Additional measurements of the same landmark only provide 
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amplifying information for an existing entry in the A matrix and thus are not added. The 
deteetion of a new landmark creates a new set of entries for the LLS process: 1) the 
position of the landmark, 2) the constraint, or line, between the AUV position and the 
landmark, and 3) the vehicle position at the time of the first measurement. The latter 
entry is defined as a pose node in community standard language. The LLS process solves 
the entire SLAM problem, consisting of all pose nodes and landmark entries since the 
mission began. This adds the benefit of being able to linearize at each pose node instead 
of trying to choose a single linearization point for all previous data as an EKF would. 

5. Mathematical Example 

To provide additional insight into the inner workings of the iSAM2 algorithm, 
consider the simplified example provided in Section IILA.3, revised slightly for clarity. 
The evolution of the iSAM2 computations will be shown over the time step to show each 
process individually. We begin by applying a coordinate frame to Figure 8, with the AUV 
starting at the origin, with a forward velocity of w = 1 m/s, no side slip velocity, and 
heading 090°, along the x-axis. For simplicity, we have set the time increment, rff, at 25 
seconds. 



Figure 8. The example problem from Section IILA.3 now with a coordinate 

system overlaid. 
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After the first time step, the AUV is now at (25,0) and has detected the feature, 
shown in Figure 9. The sonar system identifies the feature and passes the information to 
iSAM2, which runs an algorithm called Joint Compatibility Branch and Bound for data 
association to produce an estimate of the feature’s position [26]. If the feature has not 
been previously detected, iSAM2 adds the feature to its feature database and stores the 
information. Subsequent sensor measurements of the feature reduce the position 
uncertainty of both the feature and vehicle. 



Figure 9. AUV position after one 25-second time step. The AUV has detected 
the feature at [50,10] and made the pose to landmark constraint in iSAM2. 

iSAM2 solves the non-linear LLS problem to produce the estimate of both the 
AUV’s trajectory, inclusive of all pose nodes, and all features detected. Graphically, this 
can be represented through a visual sparsity pattern plot, or spy plot, which shows which 
elements of a matrix have non-zero entries. In this example, we show the A and 
R matrices to show how QR factorization supports fluid variable reordering. Figure 10 
shows the spy plot for the A matrix with the two calculated pose nodes and measurements 
of the single landmark. Figure 11 shows the R matrix following QR factorization of the A 
matrix. 
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Visual Sparsity Pattern of the A Matrix after 1 25-sec Time Step 
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Figure 10. The visual sparsity pattern of the A matrix after a single 25-second 

time step with one feature detected. 


Visual Sparsity Pattern of the R Matrix following OR Factorization 
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Figure 11. The visual sparsity pattern of the R matrix following QR factorization 
of the A matrix. Some of the entries may have changed during QR 
factorization as a result of applying Givens rotations to specific entries in the 

lower half of the A matrix. 


31 




The addition of a new measurement yields faster results when using ineremental 
updating. Instead, as previously stated, of faetoring the entire A matrix gain, the new 
measurement is simply added to the bottom of the R matrix and Givens rotations applied 
to make the new matrix, i?*, upper triangular again. In this example, that savings in 
eomputational effieieney would not be noticed, but in larger data sets, this approach 
yields significantly faster results. 

Unlike the qualitative example provided at the beginning of this section, 
recovering the uncertainty values to display them as an ellipse is not straightforward. 
Unlike an EKF, where the covariance matrix is clearly identifiable throughout the 
process, the uncertainty information is more hidden in iSAM2. In a small case such as 
this, we can recover the covariance information through from Equation (3.7). However, 
once the R matrix becomes large, this equation becomes computationally burdensome 
and a different method must be used [27]. 

P = [r^r)~' (3.7) 

In this example. Equation (3.7) produces an 8x8 matrix with the following values. 
These show the covariance information associated with the vehicle pose estimates and the 
detected landmark. It should be noted that after one time step, the uncertainty here is still 
rather small and should not be viewed as typical, especially for underwater problems 
where. 

0.01 0 0 0.01 0 0 0.01 0 

0 0.01 0 0 0.01 0 0 0.01 

0 0 0.01 -0.25 0 0.01 -0.5 0.1 

0.01 0 -0.25 6.2695 -0.0011 -0.2497 12.5105 -2.4989 

0 0.01 0 -0.0011 0.0171 0 0.0011 0.0129 

0 0 0.01 -0.2497 0 0.01 -0.5003 0.1001 

0.01 0 -0.5 12.5105 0.0011 -0.5003 25.0195 -5.0011 

0 0.01 0.1 -2.4989 0.0129 0.1001 -5.0011 1.0171 

While equation (3.8) looks cluttered, we are only concerned with, for this 
example, with the values in elements [1,1] and [2, 2], which represent the position 
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uncertainty of the AUV in both x and y. The position uneertainty of the landmark is in 
the lower right quadrant. This information is highly valuable from a mapping perspective. 

This chapter discussed the algorithmic applications of managing the stochastic 
nature of robotics. It took a close look at the iSAM2 algorithm developed by Kaess and 
Dellaert [3], [25] and illustrated the algorithm with a simple example. NPS uses iSAM2 
onboard the REMUS AUVs for researeh work and the algorithm will be the primary 
platform for integrating information from aeoustie eommunieations for collaborative 
SLAM. 
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IV. COLLABORATIVE MULTIPLE AUV SLAM 


A. RECENT WORK 

We begin the development of an MVSLAM algorithm with a review of the recent 
literature focused on this area. The existing body of work can be loosely partitioned into 
five separate areas: centralization and hierarchy, beacon-aiding, cooperative SLAM, and 
dynamic SLAM. 

I. Centralization and Hierarchy 

Early approaches to the MVSLAM problem focused on creating a centralized 
solution onboard a single vehicle, as in [28], or a joint map between all vehicles with 
common nodes that facilitates cooperative loop closure, as in [29]. Moratuwage, Vo, and 
Wang [30] focus on the creation and communication of individual submaps, which can 
then be used by other vehicles for measurements or additional data association. Finally, 
Moratuwage et al. [31] proposes a single EKE SEAM algorithm for multiple vehicles. 

Alternate approaches began to emerge as researchers postulated highly 
decentralized MVSEAM solutions, which facilitated more efficient use of processing 
power and allowed each vehicle to maintain estimates of the whole group. These 
researchers began to embrace the difficulty of having all robots communicate with all 
other robots at each of the prescribed times. Eeung, Barfoot, and Eiu [32] propose a 
framework that creates a centralized-equivalent solution in a sparsely-communicating and 
dynamic environment. Nerurkar and Roumeliotis [33] explore creating centralized- 
equivalent estimates in an asynchronous network. Bahr, Walter, and Eeonard [34] 
explored the possibility of an individual vehicle utilizing a bank of filters to track 
measurements and cooperatively localize other vehicles through trilateration. Hua et al. 

[35] present a communications-heavy approach in which all robots share their local 
sensor data. All receiving robots can then process all information and arrive at the same 
estimates, increasing robustness to individual failures. Most recently. Walls and Eustice 

[36] re frame the problem in terms of client-server relationships and transmit poses to 
allow all client vehicles to reproduce central estimates, explicitly developed for highly 
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bandwidth-limited underwater communieations. However, the authors note that their 
method does not solve for loop elosure, making the approaeh impractieal for our 
purposes. 

2. Beacon-Aided SLAM 

A classic approach to reducing position uncertainty, beacon systems can be either 
mobile or static and rely on precise position information which is supplied to robots 
through queries or constant communications. Ultra-short baseline (USBL) systems are 
very common implementations of beacon-aided navigation. This focuses on beacon- 
aiding for SLAM operations and is most closely related to the work in this thesis in terms 
of algorithmic implementation of acoustic communications. Erol et al. [37] offers an 
approach whereby mobile sensors periodically ascend to the surface for a GPS fix and 
communicate their exact position upon descent to the other nodes in the network, which 
localizes all sensors in a multi-stage algorithm. Bahr, Leonard, and Fallon [38] provide 
another very similar approach; however, the communication between vehicles and the 
follow-on algorithm serve primarily for trajectory selection vice reducing position 
uncertainty. Bahr, Leonard, and Martinoli [39] present another very similar approach that 
seeks to use one vehicle as a dedicated beacon vehicle, surfacing at proscribed intervals 
and transmitting its position information to the rest of the network. Intra-vehicle range 
estimates complete the picture. This approach is the most related to the work of this 
thesis, but is not appropriate based on the required frequency of surfacing and acoustic 
communications being counter to the principles of tactical security. 

3. Behaviors and Cooperative SLAM 

This emerging branch of MVSLAM focuses on maximizing certain information 
or other parameters through vehicle orientations, trajectories, or even path planning. 
Stipes et al. [40] offers an approach to MVSLAM that utilizes distributed control 
algorithms to adapt individual robot behaviors in real time to maximize SLAM yields and 
improve robustness. Pham and Juang [41] propose an algorithm that disperses robots 
intelligently to achieve a prescribed minimum SLAM accuracy as well as adaptively 
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balance the needs of both loealization and exploration in a oommunieations-limited 
environment. 

4. Dynamic Features in SLAM 

Perhaps the most recent branch of SLAM research, dynamic problems refer to the 
alteration of maps over time, such as in warehouse inventory problems, or the tracking of 
moving targets. Traditional SLAM architeetures only work with statie features. 
Movement of landmarks would severely degrade system performanee and potentially 
cause the operative SLAM algorithm to collapse. Lee, Clark, and Salvi [42] present a 
first-generation algorithm that ean estimate both statie and dynamie features in addition 
to the vehicle pose. They accomplish this through the use of probability hypothesis 
density filters and relate all features to the vehicle location at each time step. Abrate et al. 
[43] provides an approach that bridges traditional SLAM with the environments 
envisioned in [42] through the use of a map updating technique with the aim of long term 
mapping operations in the same physical location. 

B. PROPOSED APPROACH 

The current approaehes to MVSLAM in the literature all rely on frequent 
communications between robots to alter behaviors, estimate the full state of a group of 
robots, or transmit precise positioning information from a position source such as GPS. 
While our proposed approach will be most closely related to beacon-aided SLAM 
operations, these approaehes in the present body of SLAM literature are not appropriate 
for this work beeause of our explicit consideration of taetieal seeurity. We instead focus 
on minimum acceptable performanee, not perfeet or optimal performanee, in order to 
minimize aeoustie communieations and vehicle surfacing, thus maximizing tactical 
security. Thus, if the individual SLAM solution is aeeeptable, the AUV will not 
communicate. This heuristic ally-evaluated cost function will rely on inferential solutions 
using Bayesian methods. Additionally, unlike many of the algorithms presented in the 
previous section, we do not require the compilation of a global feature map on all 
vehicles. We instead leave the construction of the global map to post processing, where 
planners ean then decide how best to navigate through or neutralize the minefield. 
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The remainder of this seetion will proeeed as follows. First, we will derive and 
demonstrate the prineiples and utility of a Bayesian inferenee. Seeond, we will diseuss 
the value that aeoustie eommunieations ean add in redueing position uneertainty. Finally, 
we will qualitatively deseribe and develop the algorithm that will produee updated 
eovarianee estimates and show that using aeoustie communications as an additional 
measurement provides significant added value. Finally, we will discuss the performance 
metrics that will be needed to evaluate system performance. Simulation results and the 
analysis of algorithm performance will be reserved for Chapter V. 

C. INFERRING COVARIANCE 

In this section we discuss the application of a Bayesian inference as a way to 
reduce position uncertainty for n-number of AUVs. This section will present the relevant 
equations with supporting examples and provide first-order insight into the value of 
acoustic communications and ranging to uncertainty reduction, as well as explore the 
major factors that we must consider when using inferential methods. 

1. The Bayesian Inference 

Bayesian inferences use Bayes’ Theorem to update a probability estimate for a 
state as additional measurements are taken. Equation (4.1) shows Bayes’ Theorem, read 
as the probability of B given A is equal to the probability of A given B multiplied by the 
probability of A occurring all divided by the probability of B occurring. A and B can be 
any event that can be described by a probability density function. 

( 4 . 1 ) 

P{B) 

For the purposes of robotic mapping, we assume that our sensors have noise that 
can be approximated as a Gaussian or normal distribution about a given mean, p, usually 
zero, with a specified covariance, a^. The position uncertainty for an autonomous vehicle 
is initially assumed to be a Gaussian spheroid, which will iteratively change into a 
Gaussian ellipsoid as additional measurements of the surrounding environment are taken 
assuming linear dynamic and measurement models. The standard formula for a Gaussian 
distribution is as follows; 
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(4.2) 


y = f{x\/u,a)= J— g " 

^jlnor^ 

Qualitatively, a Bayesian inferenee makes an assumption about the likelihood of a 
given measurement being accurate through the use of a weighting formula that is then 
multiplied by the prior state. The following equations are evaluated using a two- 
dimensional example, which closely approximates the problem at hand. 

P{0 1 y) = ^ cP{y I e)P{e) (4.3) 

P{y) 



6 is the set of parameters that define the Gaussian distribution and y is the 
measurement. The constant, c, is a normalization constant. Since we cannot measure the 
position of the AUV directly to estimate the updated state, we are instead interested 
primarily in the inferred covariance term, which we can feed back into the host AUV’s 
iSAM2 algorithm to reduce the closed loop position uncertainty. For clarity, the variance 
of the updated distribution, denoted as the covariance at time k+1 given the covariance 
matrices at time k, is expressed as a function of n-number of input covariance matrices: 


^k+\\k+\ 


(I 

1 

/ 

\ 

-1 

( 7 ^ 

( 7 ^ 

2x2 

+ ^£1 + ., 

1 ^2x2 


= 

XX 



2:2 


/ 





(4.4) 


a. A Numerical Example 

Consider a two-dimensional problem, akin to AUV localization problems, with 
two AUVs. Their positions are not relevant to the example, but they have the following 
covariance information at time k: 
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(4.5) 


In this case, the first AUV has greater positional uncertainty than the second 
AUV, and would benefit from the second’s better localization. The first AUV needs to 
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reduce its positional uncertainty and has an acoustic communications link with the second 
AUV. The second AUV shares its covariance information at time k with the first. 
Ignoring, for the moment, the uncertain effects of acoustic communications, we can 
calculate the inferred covariance at k+1 using equation (4.4) and see how the application 
of a measurement reduces the value of the covariance matrix. Using Equation (4.4), we 
find; 
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6.6319 1.0069 
1.0069 6.6319 


(4.6) 


These numbers signify the uncertainty in the x- and y-directions as well as the 
cross-correlation between them. The result in Equation (4.6) reveals a crucial point. The 
updated result is smaller than the two prior covariance matrices, which means that any 
measurement, regardless of the precision or presence of noise, will reduce the updated 
uncertainty. 

To further emphasize the utility of the Bayesian inference, we can think of the 
covariance in a visual manner. In a Gaussian distribution, the variance represents the 
expected value of the squared deviation from the mean. The square root of the variance, 
the standard deviation, can be viewed as a confidence interval—the probability that the 
true value lies within those bounds. One standard deviation on either side of the mean 
equates to a 68.2% confidence interval, two equate to 95.4%, and three to 99.7%. Eor the 
position uncertainty of an AUV, we choose a 95% confidence interval (2a) to govern 
the size of the ellipse, meaning that the AUV has a 95% probability of being somewhere 
inside that ellipse at the given moment in time. Constructing the ellipse requires radius 
values for both the semi-major and semi-minor axes as well as the rotation angle from a 
standard Cartesian coordinate system. Erom Equation (4.4), the semi-major and semi¬ 
minor axes are the square root of the entries on the main diagonal. We find the rotation 
angle as follows [44]; 
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6'' = -tan^‘ 
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Thus, from Equations (4.4) and (4.7) we can visually represent the position 
uncertainty for our example problem in terms of 95% confidence ellipses, as shown in 
Figure 12. 

(a) Prior (b) Updated 




Figure 12. A priori 95% confidence ellipses for two AUVs centered at the origin 
in (a) and after the mathematical integration of the two through a Bayesian 
inference in (b) with the updated 95% ellipse shown in magenta. 

Figure 12 shows the benefit of this mathematical technique in terms of reducing 
the position uncertainty of an AUV. This example included two AUVs with covariance 
matrices aligned in generally the same manner. Altering the geometry to place these two 
ellipses perpendicular to each other will produce a much small updated ellipse, as we will 
see in the next section. 

2. Value of Acoustic Communications and Ranging 

This section will discuss the value of acoustic communications and ranging from 
the aspect of reducing position uncertainty. From the acoustic transmission we can gain 
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two independent measurements, with assoeiated uneertainties, to reduee the position 
uneertainty of an AUV. These measurements eonstitute the ereation of a temporary 
feature within the SLAM framework. First, the aeoustie message eontains the sending 
AUV’s state information, whieh ineludes position, heading, and position uneertainty 
(eovarianee). Seeond, the aeoustie transmission, regardless of message eontent, has a 
known time of flight, whieh we ean use to ealeulate range, as diseussed previously in 
Seetion II.C. While the WHOI mieromodem does not have the ability to diseem bearing 
from ineoming messages, we ean make a few well-founded assumptions about the nature 
of aeoustie eommunieations to exploit the additional measurement from aeoustie 
transmissions. 

For short-range navigation and traeking applieations, many different industries 
use an ultra-short baseline system, whieh eonsists of a transeeiver and transponders that 
traek vehieles using aeoustie transmissions. Most eommereial USBL systems are aeeurate 
to within 0.2% for a slant range from the transponder to the AUV and within 0.1 degrees 
in bearing, out to several kilometers. That bearing aeeuraey, already obtainable in 
multiple eommereial systems, will be the eore assumption needed to extraet the aeoustie 
transmission measurement. We will make similar assumptions with respeet to using the 
aeoustie modem in ereating a temporary feature. 

From the aeoustie message we ean aseertain the bearing of the aeoustie 
transmission by finding the range and relative bearing between the two AUVs. Using the 
ealeulated TOF and error from the DelGrosso sound speed equation, we ean determine 
the error in range as follows, where rand indieates a randomly seleeted value from a 
Gaussian distribution to aeeount for the uneertainty. 

Grange = (A.. + ’ rUUd ) T0F 

= (0-3 + 0.05 • rand) • rOF 

Assuming a bearing aeeuraey of 0.1 degrees, we ean finish eonstrueting the 
eovarianee matrix to eharaeterize the uneertainty of an aeoustie transmission: 
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a. Numerical Example Revisited 

Returning to the example from the previous section, we can show the effect of 
incorporating this second measurement. Using a baseline sound speed of 1506.16 m/s, we 
add a measure of Gaussian uncertainty from the parenthetical expression in Equation 
(4.8) , 0.343 m/s, for a final sound speed of 1506.503 m/s.3 If we assume the two AUVs 
are 1000 meters apart and bear 45° relative to each other, we calculate a TOE for the 
acoustic transmission of 0.6637 seconds. We apply the sound speed uncertainty and TOE 
to equation (4.9) to produce a covariance matrix for the acoustic measurement, in polar 
notation: 


E 


acortuns 


0.1511 0 

0 0.01 


(4.10) 


We now apply Equation (4.4) for the three matrices, the two from the previous 
section and the acoustic covariance matrix just derived, after rotating by the relative 
bearing, matrix R, to find the covariance at time k+1. Since acoustic communications 
adds additional information to the measurement (the position of the second AUV), the 
acoustic communications covariance matrix will add to the second matrix, as shown in 
Equation (4.11). Equation (4.12) shows the numerical example recalculated with acoustic 
communications. 
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3 This value was obtained by ealculating equation (2.1) with the following inputs: T = 15 °C, S = 34 
ppt, and P = 400 kg/em^. 
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(4.12) 


To better understand the magnitude of reduetion provided by the aeoustie 
transmission, we introduee the Frobenius norm, whieh provides a sealar value of the 
eovarianee matrix through a root-sum-square approaeh, and thus will show the effeet of 
noise or additional measurements more easily. 

||A|| = yJtmce(A^ A) (4.13) 


Applying this equation to the results obtained in equations (4.6) and (4.12), we 
see similarities more elearly. 
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acomms 


^+1 


9.4864 
1 = 9.5362 


(4.14) 


Equation (4.14) elearly shows that aeoustie eommunieations do not signifieantly 
alter the uneertainty information and ean thus be used to ereate temporary features in a 
SLAM framework. The example problem used in this seetion is somewhat eontrived in 
that the positional uneertainties of AUVs using SLAM algorithms will likely not be that 
large, and the governing assumptions of the aeoustie transmission must be tended 
earefully. The values were ehosen to make the results of a Bayesian inferenee explieitly 
elear. 


b. Determining Relative Performance 

The reduetion in position uneertainty with the applieation of aeoustie 
eommunieations with USBL-based assumptions is predieated on how USBL systems are 
physieally implemented. The USBL transeeiver is usually mounted on a stable surfaee 
eraft with aeeess to a preeise positioning system, sueh as GPS with an onboard INS. The 
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accurate transceiver position is what allows the relative aeeuraey of USBL solution. In 
this ease, the simulated USBL system is mounted on an AUV with not-insignifieant 
position uneertainty, meaning that the eorrelation of the aeoustie measurements to highly 
aeeurate position estimates ean no longer be assumed. However, an AUV does provide a 
stable platform from whieh to eonduet aeoustie eommunieations and ranging operations. 

The stability of the platform allows us to eonsider the relative performanee 
between vehieles as a eriterion for transmitting aeoustie messages. We formalize this 
arrangement with applieation of Kullbaek-Leibler divergenee, a eoneept in information 
theory that ealeulates the relative entropy gain between two probability density funetions 
to indieate the relative value of the two systems. Kullbaek-Leibler divergenee ean be 
expressed as shown in Equation (4.15). 


KL = 


P(x) 

q{x) 


(4.15) 


We adapt the Kullbaek-Leiber eoneept of relating probability density funetions 
examine the relative performanee between two vehieles. Equation (4.16) shows the 
relationship between two AUVs at time k. We reduee the eovarianee matrix into a unitary 
value by taking the sum of the elements on the main diagonal, or the traee. 


trace(E‘) 
trace(E‘) 


(4.16) 


Equation (4.16) does not have an explieit temporal eomponent, but it is time 
varying sinee the eovarianee matriees in both the numerator and denominator both vary 
with time as the vehieles navigate and deteet features. 


c. Aspect Dependence 

The example used thus far ineluded eovarianee matriees for the two AUVs that 
were nearly eollinear, as ean be seen visually in Eigure 12 with the semi-major and semi¬ 
minor axes of eaeh ellipse being approximately aligned. The results obtained were unique 
to the geometry of this simple problem. Consider the relative differenee between the 
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principle axis angles, whieh we ean use to determine the degree of orthogonality, A6*, 
derived from equation (4.7) and shown in equation (4.17): 


A ^ = ^2 “ ^1 ~ 


1 1 
— tan ' 
2 


2cr, 




J2 


1 

— tan 
2 


2cr„ 


^yy J 


(4.17) 


In the example we have used thus far, the evaluation of Equation (4.17) with the 
eovarianee values in Equation (4.5) yields a prineiple axis angle of 45-degrees for both 
ellipses, thus a differenee of zero, meaning the ellipses are oriented along the same 
prineiple axis. 

To show the effeet on the norm of the updated eovarianee matrix, eonsider the 
following three eovarianee matriees, with the latter two being the same size but oriented 
along opposite axes: 


21 = 
22 = 
23 = 


20 

0 

1 

0 

10 

0 


0 

1 

0 

10 

0 

1 


(4.18) 


Using Equations (4.4) and (4.13) we ean test the eombination of eovarianee 
ellipses 1-2 and 1-3 to show the full effeet of aspeet on the updated and how partieular 
geometries reduee the position uneertainty better than others. The results are shown in 
Eigure 13. 
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(a) Orthogonal Eipses 


(b) Collinear Ellipses 




Figure 13. (a) The inferential results of two orthogonal, 95% ellipses as compared 

to (b) two collinear, 95% ellipses. The updated covariance matrix in (a) 
showed significant reduction (-93% by matrix norm) whereas the reduction in 
(b) is much more modest (66% by matrix norm) and still exhibits greater 
directional uncertainty in the x-direction. 


Looking at the matrix norm indicates just how significant the reduction in position 
uncertainty is for the orthogonal case; 




20.025 


E 


1,3 

k+l\k 


= 1.3166 
= 6.6854 


(4.19) 


The orthogonal geometry produced a 93.4% reduction in matrix norm, compared 
to a 66.6% reduction in the collinear case. This simple example highlights the impact of 
geometry on position uncertainty reduction, which we can exploit to further our tactical 
security aims by only requiring aiding AUVs to transmit when the perceived reduction in 
position uncertainty for the other AUV exceeds a particular threshold. 
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D. ALGORITHM DEVELOPMENT AND INTEGRATION 


Consider the postulated operational scenario of multiple AUVs mapping a 
minefield. At present, each AUV would be operating independently of the others to 
accomplish their assigned portion of the overall mission. There is no centralized control 
of the AUVs and each AUV is responsible for minimizing its own position uncertainty, 
thus maximizing local map accuracy. In the framework proposed below, the AUVs would 
have the ability to communicate acoustically to share position uncertainty information, 
enabling all AUVs to collaboratively keep position uncertainty low for the overall map of 
the area. The basic premise considers maximizing the information gain in that if an AUV 
exceeds a specified threshold, it will broadcast a message requesting assistance. The other 
AUVs will determine if their uncertainty information can assist the broadcasting AUV 
and transmit a reply if a threshold is met. This algorithm is fully distributed and 
decentralized. This section will describe the algorithm and its construction from a broad, 
operational-level view, and the key design features that optimize the relationship between 
position uncertainty and tactical security. 

1. System Representation and Objectives 

We begin with a schematic representation. Figure 14, which shows how we will 
use acoustic transmissions to reduce position uncertainty. 
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Figure 14. A schematic representation of the proposed algorithm to reduce 
position uncertainty through acoustic communications with an emphasis on 
tactical security concerns in Steps 3 and 4. 

This recursive algorithm will operate in the following steps, external to the 
mechanics of iSAM2. This allows this algorithm to operate on a much broader set of 
SLAM algorithms. 

1. The algorithm begins with the AUV pose and information matrix 

at time k. 

2. The information matrix 7?^ is sent to this algorithm for evaluation. 

3. is converted back into a standard covariance matrix. We extract the 

position variances and calculate the Frobenius norm. The norm is then 
compared to a moving average to determine the rate of position 
uncertainty growth over time. If it exceeds a specified threshold, the AUV 
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broadcasts a message requesting assistance from any autonomous vehicles 
in the area including its position and uncertainty information. 

4. An AUV receives the broadcasted message, processes it, and determines if 
the inferential results of the two covariance matrices with the additional 
aeoustic communications measurement covariance information will 
exceed a specified threshold for relative gain. If it does, it transmits a reply 
with its position and uneertainty information. 

5. The broadcasting AUV receives the replies and processes them 
sequentially. Time domain multiple access (TDMA) considerations for 
deconflicting acoustic message transmission are assumed to be in place.^ 
The received information is processed into the correct format to infer a 
updated covariance, which is then transformed back into an information 
matrix. 

6. The updated information matrix is then reinserted into iSAM2 to be 
incorporated when the next pose node is created. 

7. The proeess repeats. 

The following subsections will elaborate on several of the design features of this 
algorithm and explain how the inferred information matrix will be integrated with 
iSAM2. 

2. Measuring Covariance Growth and the Broadcasting 
Threshold 

After the ereation of a new pose node in iSAM2, we reeover the eovariance 
matrix using Equation (3.7) and extraet the 2x2 matrix corresponding to the most reeent 
pose node ereated, which correspond to position uncertainty, to form a new 2x2 
covariance matrix. We calculate the matrix norm with Equation (4.13) and store the 
values in an array. 

At this stage we may ignore the TDMA eoneems based on the eonstruetion of the simulation model, 
but the issue of eommunieations seheduling must be examined in future work to ensure that the assumption 
holds. Changes to the TDMA strueture would not fundamentally alter the algorithm being developed here. 
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Once the length of the array ineludes thirty entries, we begin a moving average 
ealeulation between two windows of fifteen entries eaeh, using the most reeent thirty 
entries. The value was heuristieally derived. The variable t indieates the eurrent time 
step. 


window^ = mean (||Sj ||(l,t-30:t-16)^ 
window2 = mean(||Ej||(l,t -15 : 


(4.20) 


Seleeting an array length of thirty allows the initialization routines of the AUVs to 
settle out and the vehiele to travel a modest distanee before enabling the remainder of the 
algorithm. Additionally, it provides a modieum of smoothing to prevent a single 
erroneous sensor measurement from triggering an aeoustie transmission. Both aspeets 
enhanee taetieal seeurity be preventing exeessive transmissions. 

We then ealeulate the pereent ehange between the moving averages: 


AS 


lr|(t:t-30) 


window^ - window^ 
window^ 


(4.21) 


If the pereent ehange exeeeds a speeified threshold, the AUV will broadeast a 
message requesting assistanee from other autonomous vehieles in the area. We will test 
for this threshold in simulation in Seetion V. 


3. Threshold for an Informed Reply 

Using the approaeh developed in Seetion IV.C.2, all reeeiving AUVs ealeulate the 
updated eovarianee matrix using information reeeived in the broadeasted aeoustie 
message, aeoustie ranging data, and eurrent uneertainty on the host vehiele. The AUV 
then ealeulates Equation (4.16) to determine the relative gain that the broadeasting AUV 
would realize if the host AUV’s information were transmitted in reply. If the relative gain 
exeeeds a given threshold, the AUV will transmit an aeoustie message in reply. This 
approaeh reduees the number of vehieles transmitting while providing assistanee to the 
broadeasting AUV and aids in deliberate or adaptive path planning to exploit the aspeet 
dependenee properties of the Bayesian inferenee. 
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4. Acoustic Message Construction and Packet Loss 

The primer on ocean acoustics in Chapter II revealed the difficulties that we 
encounter in trying to communicate underwater. Pauli, Seto, and Leonard [45] report that 
packet loss of 20-50% is not uncommon in this domain. Given our goal for minimizing 
the number of transmissions, the acoustic messages must be constructed to be more 
robust to packet loss, thus preventing additional transmissions. This aspect of the work 
will not be explicitly tested in the simulation framework presented in Section V, but 
contributes to the discussion on the use of acoustic communications for multi-AUV 
operations. 

To account for this, each acoustic message will contain three copies of the 
information, including the timestamp, shown in Equation (4.22). Given the stochastic 
nature of acoustic communications, we heuristically deem it unlikely that two of the three 
values in any set of information will be corrupted, thus rendering the message unusable, 
otherwise the message would need to contain more than three repetitions of the same 
piece of information. While we prefer the message size to be less than 32 bytes for 
efficiency, we leave the analysis of this messaging technique to packet construction for 
future work. 


desired [^TOT ’ % ’ ] 

^^8actual ~ \_^TOT ’ ^TOT ’ ^TOT ’ % ’ % ’ % ’ ] 


(4.22) 


The receiving AUV can then process the received message and format it for 
further use in this algorithm as follows: 


msg = 


^TOA 


O (la 5 ^d'OT ^ 

Mode{Xi^,x^,xJ 


(4.23) 


The WHOI micromodem will append the time of arrival (TOA) to the message, 
along with received signal characteristics (omitted from Equation (4.23) since it is not 
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relevant to this thesis). The AUV will then take the mode of eaeh of the three information 
elements to build the message needed for follow-on proeessing. 

5. Covariance Reduction 

The estimation of the updated eovarianee proeeeds aeeording to the TDMA 
sehedule in use. We apply Equation (4.11) to eaeh message, forwarding the updated 
eovarianee after the first message into the next set of ealeulations. After eaeh applieation 
of Equation (4.11), the algorithm will store the broadeasting AETV’s pose and updated 
eovarianee information for reintegration with iSAM2. The use of iSAM2, a smoothing 
filter that eonsiders all previous data during eaeh iteration, allows for the insertion of 
these temporary features into the A matrix at the eorreet point in times, alleviating the 
framework from time lateney issues. 

6. Reintegration with iSAM2 

Eor eaeh AETV that replies, we ean treat that transmission and its assoeiated 
uneertainty as a feature node. In order to keep the map aeeurate, the eommunieations 
nodes will be maintained in a separate database from the feature nodes. This allows the 
normal NETS proeess to eontinue without requiring a new formulation of the system of 
equations. 

However, at the time of publieation, we have been unable to overeome teehnieal 
diffieulties eneountered in the iSAM2 algorithm in MATEAB despite a eoneerted effort 
to do so and extensive dialogue with Dr. Miehael Kaess. These issues stemmed from the 
eonversion of the iSAM2 eode from its native C++ environment into MATEAB that had 
heretofore gone undiseovered. As a result, we will still utilize the iSAM2 framework for 
SEAM operability, sueh as navigation and mapping. Reintegration of eovarianee 
reduetion will be demonstrated through the use of a generie eovarianee funetion that 
mirrors SEAM performanee. 


Seetion IV deseribed the algorithmie framework for using aeoustie 


eommunieations to reduee position uneertainty of a requesting AUV. It eonsidered the 
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operational level implementation and diseussed the various aspeets of the algorithm from 
a systematie approaeh. It introdueed the eoneept of a Bayesian inferenee as a way of 
fusing two measurements into an updated estimate of eovarianee and diseussed 
measuring the relative value or performanee of two eovarianee matriees using an 
applieation of Kullbaek-Leibler divergenee from information theory. 
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V. SIMULATION AND ANALYSIS 


A. OVERVIEW 

This thesis considers the value of acoustic communications in MVSLAM 
operations in reducing position uncertainty of a particular vehicle. To test the theoretical 
framework proposed in the previous section, we will conduct simulations to collect the 
data necessary to answer the following research questions. 

• How can we minimize acoustic communications while achieving 
acceptable performance in terms of position uncertainty? What constitutes 
acceptable performance in light of tactical security considerations? 

• What should the threshold be for transmitting an acoustic message, both in 
broadcast for assistance and in replying to broadcasts? 

• Can we minimize the need for GPS fixes? 

• What is the value of acoustic communications for underwater MVSLAM? 

B. EXPERIMENTAL CONSTRUCTION 

I. Operational Area 

To make this work readily adaptable for in-water testing, we have selected a 
simulation area that corresponds geographically to the marina in Monterey Bay, shown in 
Figure 15. 
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Figure 15. A satellite image of the public marina in Monterey Bay, California, 

with hypothetical vehicle tracks overlaid. 


We utilize a traditional search pattern, colloquially known as a “lawnmower” 
pattern to ensure that we fully search and map an area. We utilize four AUVs, but the 
approach described in this thesis should scale to n-number of AUVs.^ The simulation 
area shown in Figure 15 has been translated into a simulated MATLAB environment, 
shown in Figure 16. It consists of a four square kilometer area partitioned into one square 
kilometer search areas. Track spacing for the AUVs is set at 20 meters to ensure 100 
percent FLS coverage in the search area, with approximately 15 percent overlap. 


^ In theory the application of the Bayesian inference will scale to n-number of AUVs, but this will 
operationally be limited by the TDMA schedule in place. TDMA concerns are not relevant to this thesis but 
will factor into the operational implementation of this framework for large numbers of AUVs. 
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Simulation Plot for Multiple AUV SLAM with Acoustic Communications 
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Figure 16. A map of the simulated search environment with the navigation tracks 
for 4 AUVs overlaid. The red tracks indicate a search in progress. Each of the 
four labeled quadrants has a different feature density. The feature density will 

be varied across simulations. 

2. Simulation Variables 

The simulation employed in this thesis maps the notional operational area. All 
AUVs are considered to be identical to the REMUS 100 AUVs outline in Section 11.A. 

Euture work can and should consider the variation of vehicle-specific parameters, such as 
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INS or FLS accuracy, on this framework. However, for this thesis, we only consider a 
few independent variables, presented in Table 3 with brief explanations. These 
parameters will remain at the given values unless stated otherwise. 


Table 3. The independent simulation control variables as related to 
determining the performance of acoustic communications in 
MVSLAM operations. 


Variable 

Initialized Value 

Time Delay 


The time elapsed before the aeoustie 
communieations framework becomes 

aetive for the AUVs. 

200 seeonds 

Reset Time 


The time required to elapse between 
aeoustie broadeasts to allow for system 
stabilization. 

50 seeonds 

Reply Threshold 


The required reduetion in position 
uneertainty required to transmit a reply to 
an aeoustie broadcast.® 

20-90 percent 

Minimum rms Average"^ 

The minimum, average rms value 
required to enable aeoustie broadeasts. 

0.5 meters 

Moving Window Size 


The number of rms average values to be 
ineluded for pereent growth caleulations 
to trigger an aeoustie transmission. 

90 seeonds 


® The reduction is based on the application of Kullback-Leibler divergence discussed in Section 
IV.C.2.b. 

^ The rms average of the time-indexed covariance matrix norms. This parameter will be more fully 
discussed in the next section. 
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3. 


Measures of Effectiveness and Performance 


Interpreting the ehange or rate of ehange in a eovarianee matrix requires the 
eonsidered reduetion of a four-element matrix into a single value. Equation (4.13) 
provides the reduetion using the Frobenius norm. However, to monitor performanee over 
time and reduee the impaet of oseillations in algorithm performanee, we must also 
eonsider a more time-weighted metrie. To aehieve this, we ealeulate the rms average of 
the stored, time-indexed Frobenius norm values. This approaeh yields two prineipal 
benefits. First, it dampens loealized oseillations within the algorithm and thus provides a 
better indieation of system performanee over time. Seeond, it provides a meehanism by 
whieh we ean define the system performanee trade spaee and make direet eomparisons to 
INS-only solutions, with GPS fixes at various finite intervals. The rms average should be 
thought of as a single number eharaeterizing the average position uneertainty at time k. 
Going forward, the rms average of the time-indexed eovarianee matrix norms will be 
truneated in writing as the rms average. 

4. Simulation Plan 

This simulation plan ereates the framework to eolleet the empirieal data needed 
for answering the researeh questions posed in Seetion I and again at the start of this 
seetion. The answers to the researeh questions need to be derived from two distinet 
bodies of data, and thus the simulation plan is partitioned into two phases in order to test 
the appropriate variables and extraet the neeessary data. 

First, we must understand the performanee of eurrent systems, ineluding both 
SEAM and INS-only algorithms. Defining this trade spaee allows us to determine what 
eonstitutes aeeeptable performanee in this field. At a minimum, the results of this thesis 
must be no worse than the existing operational eonstruets. To establish this performanee 
baseline, we will serially vary two parameters; GPS fix interval and feature density. The 
GPS fix interval will establish the INS-only system performanee. We examine six 
different GPS fix intervals: 15, 30, 45, 60, 90, and 120 minutes. Varying the feature 
density from featureless to well-featured will establish the range of the performanee 
metries for the existing iSAM2 SEAM arehiteeture. The feature density will range from 
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0^2 features per square kilometer in 6-feature increments. The GPS data, since it does 
not depend on a random process, will be collected through two simulation runs with each 
quadrant utilizing a different fix interval. For the feature density analysis, we will 
conduct multiple simulations, varying the feature map with randomly-distributed features 
each time. Each feature density will have twelve simulation runs.^ The collection of data 
on these two parameters will define the upper and lower bounds on system performance 
needed to evaluate the contribution of this thesis. 

Second, with the results of the first phase of simulation testing in mind, we must 
systematically test the communications architecture to produce the maximum acceptable 
performance whilst maximizing tactical security by minimizing the frequency of acoustic 
communications. To accomplish this, we will vary the variables presented in Table 3 
specifically the reply threshold to produce heuristic ally-optimal performance. The 
broadcast threshold will be set on logic-based condition of two parameters; the rms 
average performance of the vehicle compared to the INS-only solution and the rms 
average performance being above the minimum rms average value. 

At this stage of the research, the qualitative nature of tactical security precludes 
the use of numerical optimization techniques and thus an analytically-derived answer for 
broadcast-reply settings will be used. The performance of the AUVs between the two sets 
of simulations will inform the tactical security analysis. These simulations will be 
conducted in environments with randomly-distributed features to significantly reduce 
unintentional correlations of system performance with a particular feature layout. 

C. SIMULATION CONSTRUCTION 

1. Performance Baseline Determination 

We begin the experimental plan by determining the trade space that AUVs 
conducting SLAM operations function in. By bounding this trade space, we allow for a 
more accurate assessment of acceptable performance for broadcast-reply thresholds. 
Table 4 shows the feature density and GPS fix interval assigned to each quadrant for each 

^ A true Monte Carlo simulation would be best, but time constraints did not permit a sizable number of 
trials to be conducted. 
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of the two simulation sets. Eaeh set will be run 12 times to ensure suffieient 
randomization of features. 


Table 4. The assignment of feature densities and GPS fix intervals for the first 
two sets of simulation runs. The runs will eolleet the performanee 
metries neeessary to define the trade spaee for SLAM operations. 


Quadrant 

Feature Density [#/km^] 

GPS Fix Interval [min] 

Setl 

1 

0 

15 

2 

6 

30 

3 

12 

45 

4 

18 

60 

Set 2 

1 

24 

90 

2 

30 

120 

3 

36 

Not Required 

4 

42 

Not Required 


To give a perspeetive on how the feature density manifests in the simulation area, 
Figure 17 and Figure 18 are provided to show the feature densities for the first and 
seeond simulation runs given in Table 4. 
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Figure 17. The simulation plot with varying feature densities by quadrant. 

Quadrant 1 is featureless, Quadrant 2 has 6 features, Quadrant 3 has 12 
features, and Quadrant 4 has 18 features. 
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Simulation Plot for Multiple AUV SLAM with Acoustic Communications 



Figure 18. The simulation plot with varying feature densities by quadrant. 

Quadrant 1 has 24 features, Quadrant 2 has 30 features, Quadrant 3 has 36 
features, and Quadrant 4 has 42 features. 

2. Threshold Determination 

The seeond aspeet of the simulation plan explores the threshold for replying to a 
broadeast for help. Our emphasis on taetieal seeurity requires this analysis. To 
aeeomplish this, we will fix the feature density and GPS intervals as shown in Table 4 
and Figure 19. 
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Table 5. The fixed simulation parameters of feature density, by quadrant, and 
GPS fix interval for the threshold determination runs. 


Quadrant 

Feature Density [#/km^] 

1 


0 

2 


12 

3 


30 

4 


42 

GPS Fix Interval: 60 min 


Simulation Plot for Multiple AUV SLAM with Acoustic Communications 



Figure 19. The simulation plot for the threshold evaluation with the given feature 
densities. Quadrant 1 is featureless. Quadrant 2 has 12 features. Quadrant 3 
has 30 features, and Quadrant 4 has 42 features. 
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This feature density arrangement was chosen to force one AUV, in quadrant 1, to 
communicate periodically. The absence of features will drive an INS-only SLAM 
solution. The remaining feature densities were selected to explore the relationship 
between feature density and the ability to effectively reduce another AUV’s position 
uncertainty. We will conduct three simulations for each of the threshold values. 

The broadcasting AUV will transmit when the average rms value is greater than 
0.50 meters and the rms average approaches or exceeds the INS-only value. As explained 
in Section IV.D.3, any AUVs receiving a broadcast for assistance will process the 
message and calculate the percentage reduction that would occur if the receiving AUV 
transmitted a reply. Therefore, we set this part of the experimental plan to evaluate the 
various transmission thresholds for an informed reply. We will evaluate in increments of 
10 percent from 20-90 percent possible reduction. 

D. RESULTS AND PERFORMANCE ANALYSIS 


1. Performance Baseline Determination 


We evaluated the GPS fix interval and feature densities using the simulation 
model described previously. This section will discuss the applicable performance 
parameters of both the GPS fix interval and feature density analyses and their 
implications for threshold determination. 

For GPS fix interval, the performance overtime will resemble a sawtooth pattern 
as the uncertainty grows over time and is eventually reduced with a GPS fix. The rms 
average of this sawtooth pattern, a partial consideration for the threshold determination in 
the next section, as a mathematically closed-form answer: 




Amplitude 


(5.1) 


We tested this analytical solution in simulation using fix intervals of 15, 30, 45, 
60, 90, and 120 minutes, shown in Figure 20. 
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Average DR RMS Values over Time for 6 GPS Fix Inten/ats 



Figure 20. The rms average of the time-indexed covariance matrix norms for six 
GPS fix intervals; 15, 30, 45, 60, 90, and 120 minutes. The values converge or 
begin to show convergence to the predicted analytical solutions for the rms 

average of a sawtooth wave. 


Figure 20 shows the rms averages converging to the analytical solution in 
Equation (5.1) as expected. This data aids in defining the trade space that we will work 
from to heuristically determine the threshold values. The GPS fix interval for the NPS 
REMUS vehicles can be manually set to any interval. We default to 30 minutes in 
practice. 

Second, we evaluated SEAM algorithm performance across eight feature densities 
to assess performance in terms of both covariance, which correlates directly with map 
accuracy, and the rms average of the covariance matrix norm, which provides a more 
stable indicator over time. Eigure 21 and Eigure 22 show the averaged results of the 12 
simulation runs in terms of covariance and the rms average of covariance, respectively. 
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Average SLAM C<rranance Main Norm [m] 


Average SLAM Covariance by Feature Den^ty for 12 Simulation Runs 



Figure 21. Average eovarianee matrix norms from 12 simulation runs in eaeh of 

eight different feature densities. 


Average Covariance Norm RMS Value by Feature Density for 12 Simulation Runs 



Figure 22. The rms average value of the eovarianee matrix norms from 12 
simulations runs in eaeh of eight different feature densities. 
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Figure 22 provides significant insight into SLAM performance across different 
feature densities. We can see that the performance does not improve appreciably once the 
number of features exceeds 30 features per square kilometer. In the future this will allow 
for adaptive search planning by multiple AUVs to ensure that search areas can be 
adaptively reassigned to allow balancing of features between AUVs. This will stabilize 
the overall position uncertainty for the AUV flight and improve map and navigational 
accuracy since each AUV will individually have better SLAM solutions following area 
reassignment. 

To determine what constitutes acceptable performance, we combined the results 
from the GPS fix interval and feature density analyses into a common plot. Figure 23 
shows the rms average of the covariance matrix norms over time for four feature densities 
and four GPS fix intervals. 


Defining the Trade Space; SLAM Covariance and INS>GPS Covariance 



Figure 23. The rms average values of the covariance matrix norms for four 
feature densities and four GPS fix intervals. The feature densities are 6, 18, 
30, and 42 square kilometers and the GPS fix intervals are 30, 60, 90, and 120 
minutes. This figure defines the performance trade space for a reply threshold 

determination. 


68 



















As stated previously, the goal is to provide aeeeptable, not optimal, performanee. 
The 30-minute fix interval represents the lower bound in Figure 23. This fix interval is 
operationally burdensome in a well-featured area, akin to the SLAM results for 30 and 42 
features per square kilometer, sinee the vehiele will have to eease mapping, surfaee, 
obtain a GPS fix, submerge, and regain traek. In the eontext of a minefield neutralization 
problem, this means the vehiele also has to reaequire the field and loealize itself to ensure 
the eorreet mine is targeted. Given that burden, we relax the GPS fix interval to 60 
minutes to provide greater flexibility and aeeommodate a wider range of feature 
densities. The rms average for the 60-minute fix interval provides the upper bound on the 
aeeeptable performanee envelope. Thus, SLAM performanee with an rms average greater 
than 0.8 meters we rejeet as unaeeeptable. 

The lower bound will govern the minimum uneertainty threshold for 
broadeasting. To prevent exeessive aeoustie eommunieations, a taetieal seeurity eoneern, 
we set the minimum rms value for transmission at an rms average of 0.5 meters or 
greater. SLAM performanee may aehieve better results than this, meaning the AUV will 
not need to eommunieate beeause the SLAM performanee alone meets or exeeeds our 
aeeeptable performanee standard. 

We will implement the upper and lower bounds algorithmieally by using logieal 
eomparisons. The broadeasting AUV will transmit an aeoustie message requesting 
assistanee from other AUVs in the area if the following two eonditions are met: 1) the 
time-indexed rms average of the broadeasting AUV’s eovarianee matrix norm is greater 
than the speeified minimum rms value of 0.5 meters; and 2) if the time-indexed rms 
average is greater than the INS-equivalent rms average minus 0.25 meters. 

An additional result to those presented in this seetion involves the standard 
praetiee of using finite GPS fix intervals. Figure 22 and Figure 23 show that, as long as 
the AUV is deteeting features periodieally, defining a modestly-featured environment, the 
performanee in terms of position uneertainty will be aeeeptable. In this ease, the GPS fix 
interval ean and should be shifted from a finite interval to a variable interval based on the 
predieted position uneertainty. This will provide greater operational effieieney by 
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minimizing the time spent ascending and descending for GPS fixes and reduce the 
probability of counter-detection when the AUV is at the surface. 

2. Broadcast-Reply Performance 

The reply threshold is based on the possible percent reduction for the broadcasting 
AUV if the replying AUV were to transmit its position and covariance information. In an 
effort to minimize acoustic communications and still achieve acceptable performance, we 
varied the reply thresholds in ten percent increments from 20-90 percent and measured 
the performance parameters. We consider the performance of the broadcast-reply 
threshold in terms of the average number of acoustic transmissions by an AUV for the 
given threshold value and by the average percent reduction achieved for the given 
threshold value. 

Intuitively, we would expect the number of transmissions to increase with a lower 
threshold for reply. However, as Figure 24 shows, the number of transmissions actually 
increases as the threshold for reply gets larger. This results from the fact that the replying 
AUVs may not be able to meet the higher threshold, leaving the broadcasting AUV 
transmitting acoustic messages more often requesting assistance. 


Average Number of Acoustic Transmissions by Threshold Value 



Figure 24. Average number of acoustic transmissions by possible percent 
reduction threshold value for each AUV. The average is across three 
simulations for each threshold value. 
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In terms of the percent reduction, we again analyze the average percent reduction 
by threshold value from the three simulations. The results are shown in Figure 25. 


Average Percent Reduction from Acoustic Communications by Threshold Value 



Figure 25. Average percent reduction from acoustic communications by possible 
percent reduction threshold value. The average is across three simulation runs 
at each threshold value. In poorly featured environments (red and blue lines), 
the trend is clear. In modest to well-featured environments (magenta line), the 
small number of simulations did not smooth the data sufficiently to draw 
conclusions. The black line represents the fourth AUV, which did not 

communicate. 


We again see overall higher percent reduction at the lower threshold, especially 
for the AUVs operating in poorly featured environments, the red and blue lines in Figure 
25. The small number of simulations did not provide sufficient smoothing of the data for 
the two AUVs operating in the well-featured environments, but their performance as 
previously been defined as acceptable and thus do not factor into this analysis. 

The decreasing trend in Figure 25 can be attributed to the greater capacity for an 
individual AUV to reply to the broadcast for assistance. As we demonstrated in Section 
IV.C.2 on the utility of the Bayesian inference, the greater number of replies will result in 
a higher overall reduction. The higher threshold value thus bars other AUVs that could 
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provide meaningful reduction from assisting the broadcasting AUV, limiting the 
effectiveness of the acoustic communications framework presented here. 

The performance of the AUVs at the 20 percent reply threshold is given in 
Figure 26. 


RMS Averages for MVSLAM with a 20 Percent Threshold 



Figure 26. RMS average values of the time-indexed covariance matrix norms for 
four AUVs using a broadcast-reply acoustic communications scheme with a 
20 percent possible reduction reply threshold. 

In comparing Figure 26 to Figure 22 and Figure 23 we can see that this acoustic 
communications framework produced acceptable performance as we previously defined 
it. AUV 4 was in a well-featured area and did not require acoustic communications to 
reduce position uncertainty. Flowever, the other three AUVs all benefitted from the 
acoustic communications framework and were able to improve their position uncertainty 
to values greater than the SLAM algorithm alone would have produced for each of the 
feature densities. 

E. VALUE OF ACOUSTIC COMMUNICATIONS 

The results of the simulations indicate unequivocally the efficacy of the acoustic 
communications framework proposed in this thesis. We recommend setting the reply 
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threshold at the lowest possible value to allow the greatest number of AUVs to render 
assistanee while minimizing aeoustic eommunieations. This is in keeping with the 
prineiples of taetical seeurity that require us to limit aeoustie eommunieations in order to 
prevent eounter deteetion of the aeoustie signals and maintain the operation eovert. 


Section V articulated the simulation model and plan and presented and discussed 
the results. We sought to answer the questions posed at the beginning of the section on 
how to minimize acoustic communications while achieving acceptable performance, 
minimizing the need for GPS fixes, and determining the informed threshold for reply in 
the broadcast-reply framework developed in Section IV. We found that acoustic 
communications do add significant value to multiple-AUV operations. Additionally, we 
explored the effect of feature density on SLAM operations and discovered the possibility 
of gaining greater operational efficiency by shifting from finite to variable GPS fix 
intervals. 


73 



THIS PAGE INTENTIONALLY LEET BLANK 


74 



VI. CONCLUSIONS AND FUTURE RESEARCH 


A. CONCLUSIONS 

1. Major Results 

The simulations and analysis of the proposed aeoustie eommunieations 
framework shows the validity of a multi-vehiele distributed system approaeh to redueing 
position uneertainty in underwater navigation. We traded optimal performanee for 
aeeeptable performanee and eonstrueted the aeoustie eommunieations broadeast-reply 
framework based on that metrie. Through simulation we determined that a low threshold 
of informed reply produees fewer aeoustie transmissions while giving the greatest pereent 
reduetion of position uneertainty. We also gained greater insight into the effeet of feature 
density on the performanee of SLAM algorithms and found a potential ability to shift 
from finite to variable GPS fix intervals. 

Additionally, the framework developed in this thesis ean be extrapolated to a 
large number of AUVs, thus providing greater flexibility to operational eommanders.^ A 
system of AUVs has the ability to share position information and eollaboratively reduee 
their own position uneertainty and inerease map aeeuraey. This inereases mission 
effeetiveness by inereasing eoverage, redueing mission time, and improving system 
robustness through the use of multiple AUVs. 

2. Contributions 

SLAM frameworks have traditionally only utilized expropioeeptive sensors for 
measurements. The use of the aeoustie modem indieates that there may other non- 
traditional sensors that ean help improve SLAM algorithm performanee by adding 
temporary features. We evaluated the SLAM algorithms through the lens of taetieal 
security. The tactical aspects of acoustic communications are an important parameter for 
multi-vehicle collaborative navigation and operations. This aspect of acoustic 
communications and SLAM has, to this author’s knowledge, never been explored before. 

^ As previously mentioned, this eonelusion is predieated on handling the TDMA seheduling eoneems. 
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We developed an information-theoretie framework eonsistent with optimal estimation for 
minimizing system navigation and mapping errors that is not reliant on external beaeon 
or positioning systems. Lastly, we developed a simulation model to address the various 
model parameters that affeet the effeetive deployment of a multi-vehiele system engaged 
in underwater navigation and mapping operations. This ineluded a brief analysis of aspeet 
dependenee as a means of aehieving greater uneertainty reduetion and is a key parameter 
for adaptive seareh and path planning in future work. 

3. Limitations and Issues 

The simple aeoustie eommunieations framework developed and presented here 
does have several limitations. First, the framework ean only sueeeed when AUVs are 
within eommunieations distanee of eaeh other. Seeond, the use of another AUV’s 
eovarianee information to reduee position uneertainty does introduee eross-eorrelation 
into the global eovarianee matrix. The introduetion of eross-eorrelation terms will 
jeopardize the independenee assumption that underpins the probabilistie framework. That 
was not aeeounted for in this thesis and must be analyzed as a matter of future work. We 
did not eonstruet the arehiteeture to allow for other reply AUVs to reduee their position 
uneertainty from other reeeived replies. The omnidireetional nature of aeoustie 
eommunieations means that eaeh AUV that replies has the ability to reeeive the other 
replies and use those messages to reduee their own position uneertainty, even if it meets 
or exeeeds the aeeeptable performanee standard defined in this thesis. 

B. FUTURE WORK 

This thesis eonsidered a very narrow segment in the nexus between aeoustie 
eommunieations and underwater MVSLAM operations. Given that narrow seope, we 
neeessarily leave several next steps to future researehers. 

• How is the independenee of eovarianee matriees affeeted by aeoustie 
eommunieations between vehieles? 

• How ean we use this implementation to induee loop elosure in SLAM? 
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• What are the proper contents for the acoustic message to keep the size 
under 30 bytes while still achieving the required single-transmission 
success metric? 

• How can this implementation evolve adaptive path planning and mapping? 

• How can we leverage any transmissions to create a communications map 
of the environment? 

• How can we use acoustic communications to distribute highly accurate 
position information to a group of AUVs, such as when a single AUV 
obtains a GPS fix? 

C. APPLICATIONS 

The algorithm for acoustic communications in support of underwater MVSLAM 
operations applies specifically to the operational environments detailed in Section I, but 
we can also draw greater value by expanding the view of this thesis to a much broader 
level. This thesis implemented the novel consideration of an acoustic transmission as an 
independent measurement, not simply a communications path. The value taken from both 
the transmission itself and the information contained in the transmission enabled us to 
reduce the position uncertainty of an AUV, thus enabling greater operational flexibility 
and enhancing tactical security. Taking this idea more abstractly, we can see the inherent 
value of using non-traditional sensors to better localize a ship’s position in an A2/AD 
environment. In the future, this may consist of a ship launching an AUV to obtain a 
precise position fix in a well-featured or contoured environment that the ship cannot enter 
into or a submarine utilizing other sensors, such as a fathometer, mapping sonar, or 
periscope, to obtain position information and thus reduce position uncertainty. The ability 
of our forces to sufficiently localize their position in an A2/AD environment will be a 
crucial warfighting requirement in the future, and conceptualizing various non-traditional 
sensors as possible sources of position information constitutes a modest first step in that 
direction. 
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APPENDIX A: SOUND SPEED EQUATION CONSTANTS 


C„ 0.50110939883xl0‘ 

Q 2 -0.550946843172x10' 

C„ 0.221535969240x10^' 

C 51 0.132952290781x10' 

C ,2 0.128955756844x10 ' 

Cp, 0.156059257041x10° 

Cp 2 0.244998688441x10 " 

Cp, -0.883392332513x10 * 

Cps -0.127562783426x10' 

Cpp 0.635191613389x10 ' 

Cp 2 P 2 0.265484716608x10 ' 

Cpp2 -0.159349479045x10 ' 

Cpp, 0.522116437235x10 ° 

Cp,p -0.438031096213x10 ° 

Cs2P2 -0.161674495909x10 * 

Cp 2 s 0.968403156410x10^ 

Cps 2 P 0.485639620015x10 ' 

Cpsp -0.340597039004x10 ' 
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